###########读取文件并进行INLA分析

######ASE function要求结果file文件为count数，且顺序为
######twin1_control_ref,twin1_case_ref,twin2_control_ref,twin2_case_ref,
######twin1_control_var,twin1_case_var,twin2_control_var,twin2_case_var.

library(data.table)
library(magrittr)
library(parallel)
library(dplyr)
library(boot)
library(INLA)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/varscan")
sel1=c("M30_M29","M26_M25","M35_M36","M28_M27")
sel1=c("X2B_X1T","M8_M7","M6_M5")
sel1=c("M2_M1","M48_M47","M50_M49")
sel1=c("M18_M17","M20_M19","M22_M21","M40_M39")
tag=1
sel1=paste0(sel1,".snp")


ASE=function(x)
{
if (length(which(!is.na(x)))/4>=1){
sel=which(!is.na(x))
rct=x[sel]
ref=rct[1:(length(rct)/2)]
var=rct[(length(rct)/2+1):(length(rct))]
if (length(sel)/4>1){
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4),x2=rep(c(0:(length(sel)/4-1)),each=2))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1","x2")
#colnames(df)=header
formula = y ~ 1 + f(x2, model = "iid") + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1 + f(x2, model = "iid")
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}else{
df=data.frame(y=var,Ntrials=ref+var,x1=rep(c(0,1),length(sel)/4))
#header=c(rep("y",length(sel)/2),rep("Ntrials",length(sel)/2),"x1")
#colnames(df)=header
formula = y ~ 1 + x1
m1 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
formula = y ~ 1
m0 <- inla(formula, data = df, family = "binomial", Ntrials = Ntrials,quantile = c(0.005, 0.025, 0.975, 0.995))
b10 <- exp(m1$mlik[2] - m0$mlik[2])
}
return(b10)
}
}

for(i in 1:length(sel1)){
file=read.table(sel1[i],head=T,sep="\t")
file=file[file$exits=="T",]
ref=paste(c("normal","tumor"),"_reads1",sep="")
var=paste(c("normal","tumor"),"_reads2",sep="")
file1=file[,c(ref,var)]

rt=apply(file1,1,ASE)
result=as.numeric(as.character(rt))
file$BayesFactor=result
file_name=paste0("../bayes_BF/",gsub(".snp","",sel1[i]),".bayes_factor.txt")
write.table(file,file_name,quote=F,row.names=F,sep="\t")
}